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We present a detailed first principles study of Fe-pnictides with particular emphasis on competing 
magnetic interactions, structural phase transition, giant magneto-elastic coupling and its efltect on 
phonons. The exchange interactions Jij(R) are calculated up to ~ 12 A from two different 
approaches based on direct spin-flip and infinitesimal spin-rotation. We find that Jij{R) has an 
oscillatory character with an envelop decaying as 1 /R^ along the stripe-direction while it is very short 
range along the diagonal direction and antiferromagnetic. A brief discussion of the neutron scattering 
determination of these exchange constants from a single crystal sample with orthorhombic twinning 
is given. The lattice parameter dependence of the exchange constants, dJij /da are calculated for 
a simple spin-Peierls like model to explain the fine details of the tetragonal-orthorhombic phase 
transition. We then discuss giant magneto-elastic effects in these systems. We show that when 
the Fe-spin is turned off the optimized c- values are shorter than experimetnal values by 1.4 A 
for CaFe2As2, by 0.4 A for BaFe2As2, and by 0.13 A for LaOFeAs. We explain this strange 
behavior by unraveling surprisingly strong interactions between arsenic ions, the strength of which 
is controlled by the Fe-spin state through Fe-As hybridization. Reducing the Fe-magnetic moment, 
weakens the Fe-As bonding, and in turn, increases As-As interactions, causing a giant reduction 
in the c-axis. These findings also explain why the Fe-moment is so tightly coupled to the As-z 
position. Finally, we show that Fe-spin is also required to obtain the right phonon energies, in 
particular As c-polarized and Fe-Fe in-plane modes that have been recently observed by inelastic 
x-ray and neutron scattering but cannot be explained based on non-magnetic phonon calculations. 
Since treating iron as magnetic ion always gives much better results than non-magnetic ones and 
since there is no large c-axis reduction during the normal to superconducting phase transition, the 
iron magnetic moment should be present in Fe-pnictides at all times. We discuss the implications 
of our results on the mechanism of superconductivity in these fascinating Fe-pnictide systems. 

PACS numbers: 74.25. Jb,67.30.hj,75.30.Fv,75.25.tz, 74.25. Kc 



I. INTRODUCTION 

The recent discovery of superconductivity at Tc's up to 
55 K in iron-pnictide systemsii^ii^iii has sparked enormous 
interest in this class of materials. So far four types of 
materials have been discovered. The first one is the rare- 
earth pnictide oxide layered systems, REOFeAs which 
is denoted as " The second class is the so 
called " 122" systems with the chemical formula MFe2As2 
(M=Ca,Sr, etc)^^^^^^. The third system with = 18 K 
is MFeAs (M=Li and Na), which is similar to REOFeAs 
but instead of REO-layers, we have now small alkali met- 
als such as Li—. The last one is the binary Fe(Se,Te) 
systems which have been shown to superconduct up to 
12 K under pressureii. 

The crystal structures of these four systems are shown 
in Fig. 1. The common feature in these Fe-pnictide su- 
perconductors is the presence of FeAs plane (or Fe(Te,Se) 
in the case of 11 systems), which is shown in Fig. 2. Ba- 
sically Fe atoms form a regular square lattice just like 
the Cu02 plane in cuprates. However the important dif- 
ference is the location of the arsenic ions which are not 
located between two Fe ions but rather above/below the 
center of Fe-square. This arrangement of arsenic ions 
has several important consequences in the electronic and 



(a) 11 11 (LaOFeAs) (b) 122 (CaFe^As,) (c) 1 11 (LiFeAs) 




FIG. 1: (color online) The crystal structures (with origin 
choice 1) of four types of Fe-pnictide systems that have been 
discovered so far. 



magnetic properties of these systems. Since As is not di- 
rectly between two Fe ions, the Fe-Fe distance is not large 
and direct Fe-Fe overlap plays an important role in the 
band formation near the Fermi level. Then, the delicate 
interplay between Fe-Fe, Fe-As, and even As-As inter- 
actions (which is very important in 122 systems such as 
Cal22) result interesting electronic and magnetic prop- 
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erties that are super-sensitive to the As-z position and 
the c-lattice parameter of the Fe-pnictide system. In this 
paper, we wiU focus on the structural, dynamical and 
magnetic properties of 1111 and 122 systems only. For a 
recent review of electronic and superconducting proper- 
ties of Fe-based superconductors, we refer the reader to 
David Singh'si^^ and Igor Mazin's^- articles in this issue 
and references therein. 

A common phase diagram for iron-pnictides has 
emerged^* in which the stoichiometric parent compound 
shows a structural anomaly around 150-200 K, be- 
low which spin-density-wave (SDW) antiferromagnetic 
ordering^i^ i^^i^^ appears, which is due to nesting Fermi 
surfaceaiiii^ii^i^S. The SDW ordering is further stabi- 
lized against the normal checkerboard antiferromagnetic 
ordering (denoted as AFl) due to strong antiferromag- 
netic interactions along the Fe-square diagonai^i. Su- 
perconductivity in these systems only occurs when the 
SDW ordering and the structural distortion are sup- 
pressed, which can be achieved in a number of ways 
such as fluorine doping on the oxygen siteii>^>^i^, or 
hole doping (Lai-xSix)^'^^^'^ or by applying external 



The structural distortion which is common to all par- 
ent compounds can be characterized by either primitive 
monoclinic space group P2/c (P112/n) or the conven- 
tional orthorhombic cell with space group Cmma^^.. The 
relation of these two representations is indicated in the 
bottom panel of Fig. 2. We note that when the sys- 
tem is distorted (i.e. 7 ^ 90.0) one of the Fe-pairs gets 
closer and the other Fe-Fe distance gets longer, yielding 
an orthorhombic lattice (i..e Qq 7^ bo)- Below we will 
successfully explain how this structural phase transition 
is tightly coupled to the magnetic SDW ordering. Fi- 
nally we note that the reported space groups in the SDW 
statei^!?^ (i.e. Cmma or P2/c) are actually the space 
groups of the system without the Fe-magnetic moment. 
Hence technically the Cmma is not the right space group 
for the SDW magnetic system. If we ignore the antifer- 
romagnetic stacking of the FeAs planes, the actual space 
group is Pbnm which is primitive as expected. This will 
be important in the discussion of magnetic phonon calcu- 
lations in Sec. VI where we can not use the space group 
Cmma but has to use Pbmb for 1111 systems. The sit- 
uation is very similar for the 122, 111, and 11 systems 
as well. For example, the Fmmm space group of Ill- 
systems is actually reduced to Bbmb (spg. number=66, 
origin 1) when the iron-spins arc considered. 

Clearly, the understanding of electronic, magnetic, and 
structural properties of the parent FeAs compound is the 
key to determining the underlying mechanism that makes 
these materials superconduct upon electron/hole doping. 
In this paper we present a detailed first-principles study 
of Fe-pnictides with main focus on the competing mag- 
netic spin-interactions, structural phase transition, the 
giant magneto-elastic coupling and the phonons. Our 
main objective is to demonstrate that Fe-spin is the key in 
understanding many properties of these systems, includ- 




FIG. 2: (color online) Top: A view along c-axis of the FeAs- 
plane and the relations between the primitive and \/2 x \/2 
supercell used in our calculations. The dark and light shaded 
areas indicate the As atoms below and above the Fe-square 
lattice, respectively. Bottom: Relation between conventional 
(Cmma) and primitive (P2/c) cells of the orthorhombic struc- 
ture. 



ing lattice parameters, atomic positions, and the phonon 
spectrum. When the Fe-spin is ignored and non-magnetic 
calculations are done, the results do not agree with most 
of the experimental data. This observation could be the 
key in identifying the mechanism of superconductivity in 
these systems. 

This paper is organized as follows. In the next section, 
we discuss the energetics of possible spin configurations in 
Fe-pnictides within a unified model from all-electron fix- 
spin moment calculations. We will show that the SDW 
magnetic ordering is the only stable ground state for Fe- 
pnictide. In Sec. Ill, we will calculate the exchange in- 
teractions Jij{R) up to ~ 12 A using two different ap- 
proaches based on direct spin-flip and infinitesimal spin- 
rotation. We find that Jij (i?) has an oscillatory charac- 
ter with an envelop decaying as l/R^ along the stripe- 
directions. On the other hand, it is short range along 
the diagonal direction and antiferromagnetic, suggesting 
it is superexchange type and an important contributer 
towards the stabilization of SDW ordering. A brief dis- 
cussion of the experimental determination of these ex- 
change constants from an orthorhombic-twin crystal is 
also given in this section. In Sec. IV, we will discuss 
the tetragonal-orthorhombic lattice distortion. We will 
calculate the lattice parameter dependence of the ex- 
change constants, dJij/da, and then use it in a simple 



3 



spin-Peierls like model to explain the fine details of the 
tetragonal-orthorhombic phase transition that is driven 
by the SDW ordering. In Sec. V, we will discuss the gi- 
ant magneto-elastic effects in these systems where iron- 
spin controls the strength of Fe-As and As- As hybridiza- 
tion which results huge dependence of the magnetic and 
structural properties on the As-z and c-axis of the lat- 
tice. Finally, in Sec. VI, we show that Fe-spin is also re- 
quired to obtain the right phonon energies, in particular 
As c-polarized and Fe-Fe in-plane modes that have been 
recently observed by inelastic x-ray and neutron scatter- 
ing measurements but could not been explained based on 
non-magnetic phonon calculations. Our conclusions will 
be given in Sec. VII. 
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II. SPIN DENSITY WAVE (SDW) ORDERING 

The early theoretical studies identified several can- 
didate ground states for Fe-pnictides such as a non- 
magnetic metal near a ferromagnetic or antiferromag- 
netic instabilityiii^i^ and a simple antiferromagnetic 
semi-metal22ii2^. The later calculations suggested that 
Fe-pnictide has an antiferromagnetic spin-density-wave 
(SDW) ground state^^ that is stabilized by Fermi-surface 
nestingi2ii2£ as well as by strong antiferromagnetic spin 
interactions along the Fe-square diagonal^^. Recently 
several nice review articles have been published^i^ 
about the Fermi surface nesting, band structure prop- 
erties, and the delicate interplay between the structural 
and electronic properties'^ in these fascinating systems. 
Hence, here we will not discuss the Fermi-surface nest- 
ing and band structure properties but rather focus on 
the energetics of different spin-configurations in order to 
determine the nature of magnetic interactions present in 
these systems. 

In order to demonstrate the SDW ground state of Fe- 
pnictide systems, one needs to consider V2x \/2 supercell 
of the tetragonal cell shown in Fig. 2 with four different 
magnetic spin configurations. These are non-magnetic 
(NM, i.e. no spin polarization), ferromagnetic (F) and 
the two different antiferromagnetic spin configurations 
shown in Fig. 3. The first one of the antiferromagnetic 
configurations is AFl where the nearest neighbor spins 
are anti-parallel to each other. The second antiferromag- 
netic configuration, AF2, is shown in Fig. 3c and 3d. In 
AF2 the Fe spins along the square diagonal are aligned 
antiferromagnetically. This is the stripe-phase which was 
first predicted from Fermi-surface nestin g^^'^" . The AF2 
spin configuration can be considered as two interpenerat- 
ing simple square AF sublattices (circle and square sub- 
lattices in Fig. 3c). From the classical Heisenberg ener- 
gies of AFl and AF2, one sees that the AF2 spin config- 
uration is stabilized when J2 > Ji/2. We note that in 
AF2 spin configuration, due to large antiferromagnetic 
J2 interactions, the spins along the diagonal direction arc 
aligned antiferromagnetically, forcing spins to be parallel 
and antiparallel along the a— and 6— directions. Hence 




FIG. 3: (color online) Four possible magnetic configurations 
for the Fe-square in Fe-pnictide and the corresponding energy 
expressions in terms of a simple J1-J1-J2 model. Two antifer- 
romagnetic configurations are considered in this study. Top- 
right panel (b) shows the AFl configuration where nearest 
neighbor spins are always aligned anti-parallel. Two bottom 
panels (c-d) show the AF2 configuration where the next near- 
est neighbor spins (i.e., J2) are always aligned anti-parallel. 
Note that this is the same stripe-phase predicted from Fermi- 
surface nestingi^i^ and it is frustrated. 



the Ji interaction can not be fully satisfied and therefore 
the system is called frustrated. In frustrated magnetic 
systems, it is known that the frustration is usually re- 
moved by either a structural distortion or by other ef- 
fective spin-spin interactions that originate from thermal 
and quantum fluctuations of the spins^i^^. As we shall 
see below, in Fe-pnictide the Heisenberg picture is only 
an approximate model. The exchange interactions de- 
pend on the spin-configurations considered and therefore 
Ji could be even ferromagnetic in AF2, removing the 
frustration totally. Finally it has been recently shown 
that the total energy of the AF2 spin configuration in- 
creases as sm{9)^ when the two interpenetrating AFM 
sublattices are rotated by an angle 9 with respect to each 
other'^^. Hence the infinite degeneracy of the classical 
Heisenberg model of AF2 structure has been already re- 
moved by either non-Heisenberg like interactions and/or 
long range interactions that are present in Fe-pnictides. 

In order to determine which spin configurations among 
NM, F, AFl, and AF2, is the ground state, we have 
carried out total energy calculations for each case us- 
ing experimental structure. The calculations were done 
using the full-potential linearized augmented plane- wave 
(FP-LAPW) method, within local density approximation 
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FIG. 4: (color online) (a) The total energy per Fe atom versus 
magnetic moment for F, AFl and AF2 spin-configuration, 
indicating AF2 is the only ground state of the system, (b) 
The magnetic interactions for nn and nnn Fe ions obtained 
from the energies of F, AF, and AF2 configurations. 



(LDA) using Perdew-Wang/Ceperlye-Alder exchange- 
correlation^Sdi,. We also used the ultrasoft pseudo po- 
tential planewave (PW) method^ for cross checking of 
our results and for phonon calculations. Since in spin- 
polarized calculations it is very easy to get a local min- 
imum, we followed a different strategy. In our calcula- 
tions we fixed the magnetic moment per Fe ion and then 
scanned the total energy as a function of Fe-magnetic 
moment. Our results are summarized in Fig. 4. The 
zero of energy is taken as the M=0 case (i.e., NM calcu- 
lation). From Fig. 4, it is clear that LaOFeAs has only 
one magnetic ground state which is AF2. The Ferro spin- 
configuration always results the highest energy regardless 
the Fe-magnetic moment. Similarly AFl ordering always 
yields energies higher than the NM case. For the AF2 or- 
dering, we see that the energy minimum occurs near the 
fixed moment calculation with M=l. Repeating calcula- 
tions where magnetization is not fixed, we obtained the 
optimum magnetic moment as M=0.87 [Ib per Fe. As 
we discuss below in detail, the Fe magnetic moment is 
further reduced almost by a half when the structure is 
allowed to distort to due to AF2 stripe ordering. 

One confusion with the DFT studies of Fe-pnictide is 
the calculated Fe-magnetic moment. The most of the 
calculations, in particular those based on pseudopoten- 
tials, give too large moment around 2.0 hb compared 
to experimental values of 0.3 — 0.8 \ib- In order to ex- 
plain the small experimental moment, several theories 
based on quantum or thermal fluctuations have been pro- 
posed since the SDW system is magnetically frustrated'*^. 
On the other hand, there are studies suggesting that the 
small moment is due to electronic effects, local chemistry 
of Fe and its interaction with the As ions. It has been 
shown that a small displacement of As-z position and/or 
structural distortion can easily change the Fe-moment 
from 2.0 \ib to 0.5 /x^i^ii^iii. In section V, we will 
discuss this high sensitivity of the Fe-moment to the As- 
z position in detail in the context of giant magneto-elastic 
couping since z-position and c-axis of the crystal are cou- 
pled and have to be treated equally. 

In order to gain a better insight into the nature of the 



magnetic interactions present in Fe-square lattice of the 
Fe-pnictide system, we map the calculated total energies 
of the F, AFl and AF2 configurations shown in Fig. 4a to 
a simple Heisenberg like model H = Ep -I- ^ JijMiMj 
for a given fixed Fe moment Mi . For fully localized spin- 
systems this is a perfect thing to do but for the case of Fe- 
pnictide this is only an approximation. Nevertheless, the 
calculated Js should be a good indication of the magnetic 
interactions present in the system. We also note that 
these interactions are valid at high temperatures above 
the magnetic ordering transition where spin-flips are the 
relevant magnetic excitations (and not the spin- waves). 
Fig. 4b shows the effective Ji and J2 obtained from the 
energies of the F, AFl and AF2 at given magnetic mo- 
ment. The dependence of the Js on the magnetic moment 
further suggests that the simple Heisenberg Hamiltonian 
is not a good model for this system. Here we use the 
term effective because the calculated Ji j is actually the 
interaction between spins i and j plus the infinite sum of 
interactions between their periodic images. The calcula- 
tions of magnetic interactions up to 4r(i-nearest neighbor 
will be discussed in the next section. From Fig. 4, it is 
clear that both Ji and J2 are quite large and positive 
(i.e., antiferromagnetic). J2 is always larger than Ji/2 
and therefore AF2 structure is the only ground state for 
any given moment of the Fe ion. By looking at the ex- 
change paths for Ji and J2 (shown in insets to Fig. 4) , we 
notice that the Fe-As-Fe angle is around 75° and 120° for 
nn and nnn Fe-pairs, respectively. Hence it makes sense 
that the 2nd nn exchange interaction is as strong as the 
nn exchange because the angle is closer to the optimum 
value of 180°. 

We note that there are now a large number of studies 
of exchange interactions in Fc-pnictides based on vary 
different methods such as mapping energies to a Heisen- 
berg model^ii^ii^, linear response theories^^^i^, strong 
coupling perturbation calculations of superexchange in- 
teractions within a tight binding model^ and Kugel- 
Khomskii type effective Hamiltonian with spin and or- 
bital degrees of freedom'*^.. All these studies indicate 
that the major exchange interactions in Fe-pnictides are 
Ji and J2 that are large, comparable in magnitude and 
antiferromagnetic (i.e. frustrated). This seems to be the 
intrinsic property of FeAs plane that is common to all 
Fe-pnictide superconductors. It is quite surprising and 
also very interesting that there are strong and competing 
antiferromagnetic interactions in the Fe-pnictide system 
that result in a totally frustrated AF2 spin configuration. 
This is very similar to the magnetic ground state of the 
cuprates where the AF ordered 2D square lattices of the 
adjacent planes are frustrated^!. Even though electron 
doping seems to destroy the long-range magnetic order, 
the short range spin fluctuations will be always present 
and probably play an important role in the supercon- 
ducting phase, much like the high Tc cuprates. 

Finally we note that there have been several neu- 
tron scattering measurement o^^i^^'^^i^^ i^ of the spin- 
wave spectrum in 122 systems indicating Ji -I- 2J2 ~ 
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100 ± 20 meV. From Fig. 4, we get Ji + 2 « 80 meV for 
S—1 in 122 systems. Noting that our calculations were 
done much before the experimental measurements, the 
agreement is quite good and give confidence that DFT 
calculations actually work for predicting properties of Fe- 
pnictides. Similarly, based on our calculations^, we had 
also predicted that orthorhombic lattice parameter along 
the parallel-aligned spin-direction should be shorter than 
the axis along the anti-parallel spin direction in the SDW 
phase, which has now been confirmed by experiments^^. 
This further assures that all-electron DFT calculations 
capture many fine details of the physics in Fc-pnictides. 



III. COMPETING MAGNETIC INTERACTIONS 
IN FE-PNICTIDES 

In previous section we estimated some effective ex- 
change interactions from the total energies of F, AFl, 
AF2 but from those estimates it is not clear if the cal- 
culated parameters are limited to the nearest neighbor 
interactions or not. If the exchange interactions origi- 
nate from Fermi-surface nesting then they should be long 
range. If they originate from Fe-As-Fe superexchange 
then they should be short range. Hence by calculating 
Jij {R) as a function of Fe-Fe distance R, we can deter- 
mine if the Fermi surface nesting is the major factor in 
the magnetic exchange interactions and learn more about 
nature of these interactions and the way they couple to 
the lattice and the structural phase transition. 

Here we address this issue by calculating exchange pa- 
rameters from two different methods. We developed a 
systematic approach where the exchange parameter be- 
tween spin-« and spin-j is obtained from the total energies 
of a reference magnetic configuration and those configu- 
rations obtained by flipping the spins i and j one at a 
time and simultaneous flipping both spins. From these 
four energies, it is possible to obtain the exchange con- 
stant between spin i and j. The details of this technique 
are presented in Appendix A. From now on, we refer 
to this method as "direct spin-flip method". Note that 
this method is more appropriate at very high tempera- 
tures where the relevant spin-excitations are spin-flipping 
rather than spin- waves. 

The second method is based on linear response 
perturbation theory using Green function approach 
within rigid-spin-approximatio n^^i^^i^^ as implemented 
in openmx package^. Basically, we calculate the small 
energy change due to an infinitesimal rotation of spin i 
and spin j using Green function perturbation theory with 
the assumption that the magnitude of spins are fixed and 
only their orientations are changed. This is a question- 
able approximation for Fe-pnictides because the Fe-spin 
magnitude is very sensitive to the spin pattern as well as 
the As-position. In Sec. II, we have already seen that 
the spin-magnitude goes to zero if they are forced to be 
aligned ferromagnetically. With this in mind, we still 
think that when the system is at very low temperature 




FIG. 5: (color online) Schematic representation of a spin- 
interaction energy versus the angle between spins i and j for 
a Heisenberg (black) and non-Heisenberg model (red). The 
linear response theory gives the second derivative of the en- 
ergy at the local spin-configuration while the spin-fiip method 
defines the exchange constant as the energy difference between 
parallel and antiparallel spin-configurations. For the Heisen- 
berg model (black) both methods give the same energy, i.e. 
J. However for non-Heisenberg model, the linear response 
theory gives K (i.e. ferromagnetic for this example) while 
spin-flip method gives J (i.e. antiferromagnetic) . Needless 
to say, both are correct! K describes the spin-dynamics near 
the bottom of the spin-interaction potential (i.e. spin-waves) 
while J describes the spin-dynamics near the energy barrier 
(i.e spins are flipping as in the paramagnetic phase). 



where the spin-moment is fixed and spin- waves are valid, 
this perturbation approach should give physically correct 
results. Similar approach have been successfully used to 
study exchange interactions and spin-wave stiffness con- 
stants in transition metals such as bcc Fe^. 

It is very important to note that the calculated quan- 
tities from spin-flip and linear-response-theory are not 
the same thing and therefore these two methods can 
give totally different results for non-Heisenberg spin- 
interactions. This is schematically shown in Fig. 5. The 
linear response-theory gives the second derivative of the 
spin-interactions from infinitesimal rotation of the spins 
near their equilibrium positions. On the other hand, 
the spin-fiip method gives the energy difference between 
parallel and antiparallel spin configurations. When we 
deal with Heisenberg interaction, J{9) is proportional 
to cos(6') and therefore both the second derivative and 
the energy difference give the same answer. However 
for a non-Heisenberg model where the dependence of J 
on the angle is not a simple cosine as shown in Fig. 5, 
the linear response theory will give the 2nd derivative 
at the local structure, i.e. K in Fig. 5 while the spin- 
flip method will give J and therefore the two methods 
will differ. This should not be considered as one of the 
method is not working but rather as an evidence that 
the spin-spin interaction is not a traditional Heisenberg 
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FIG. 6: (color online) The 3^(2) x3\/(2) x 1 supercell consid- 
ered in calculations of the exchange parameters in real-space. 
For simplicity, we show only the Fe ions with their number 
labels and spins (+ up, - down) according to AF2 ordering. 
The 3i,j for each iron pairs i and j are plotted in Fig. 6. The 
shaded region indicates the half of the supercell in which we 
can extract the exchange interactions between any pairs of 
spins (which are numbered from 1 to 16 for convenience). 



type. As we shall see below, this is indeed the case for 
Fe-pnictide where for the parallel-spin direction in the 
stripe phase, linear response theory will give a ferromag- 
netic interaction while the spin-flip method will give an 
antiferromagnetic interaction. Since the linear response 
theory probes the energy changes near the local-spin 
configurations due to small spin-rotations, the Js from 
this method should be used to explain the experimental 
spin-wave spectrum. However for spin-dynamics near or 
above the paramagnetic phase transition the spin-waves 
are not valid and therefore the Js from spin-flip method 
is more appropriate. Of course, since the Js depend on 
the spin-configurations considered one can not study the 
phase diagram of Fe-pnictides from zero to high tempera- 
tures with a single Heisenberg like Hamiltonian. In that 
case, it is important that one goes beyond the Heisen- 
berg picture and consider more complete models such as 
Kugel-Khomskii Hamiltonian where the orbital and spin- 
degrees of freedom are treated self consistently on equal 
footing^ . 

Finally, we note that both methods discussed above 
are based on periodic supercell approach and therefore 
what we calculate is actually the sum of the interactions 
between spin i and j and the interaction between their 
periodic images. Hence it is important to consider a very 
large supercell to make sure what we get is actually the 
individual exchange constant between spin i and spin j. 
Hence, we consider 3-\/(2) x 3y^2) x 1 cell which has 
total 144 atoms, 36 of which are iron. Fig. 6 shows the 
labeling of the 32 iron atoms in the large supercell. The 
direct spin-flip calculations were done using the plane- 
wave code pwsci^ with cutoff energy of 30 Ry and charge 
cutoff of 240 Ry using PBE-GGA exchange functional. 
We used 2x2x3 k-points. The perturbation calculations 
with rigid-spin-approximation were carried out using the 
package openmx^^ which implements numerical atomic 



arbitals and norm-conserving pseudopotentials. We used 
equivalent cutoff, k-point grid and the same PBE ex- 
change functional as the pwscf calculations. We used 
experimental atomic positions for the LaOFeAs system 
without any structural relaxation. Both methods (i.e. 
pwscf and openmx) give an iron magnetic moment close 
to ~ 2^B , a typical value obtained from pseudo-potential 
based methods. For convenience, in the discussion below, 
we report JS^ by setting S=l (rather than taking S=2 
and recalculating Js). 

Our results, for the exchange constants, Jij{R) from 
both spin-flip and perturbation methods are shown in 
Fig. 7. We first discuss top panel in Fig. 7 which is from 
direct spin-flip method. The strongest interactions, i.e. 
Ji.2,Ji,3, and Ji_4 are all antiferromagnetic and compa- 
rable to each other, consistent with our previous results 
obtained from the total energies of three spin configura- 
tions discussed in section H. We note that the symmetry 
between a and b axis is broken and we obtain different 
(i.e. Ji,2) and (i.e. Ji.s). However the difference 
is small and they are both antiferromagnetic. The J2I Ji 
is around 1 (here J2 is J1.4) and therefore the AF2 (i.e. 
SDW) is the ground state. Fig. 7 also shows how the 
exchange constants decay with distance along the Fe-Fe 
square diagonal, and along a- and b-directions. We note 
that the anti-parallel-spin alignment is taken to be along 
the a-axis. Interestingly, Ji j [R) already changes sign and 
become ferromagnetic at the 2nd (i.e. Ji^s) and 3rd shells 
(i.e. Ji,io) along the o— and 6-axes, respectively. As 
shown in Fig. 7, along the a— and 5— directions, Jij{R) 
has an oscillatory character with an envelop decaying as 
1/R'^. One exception to this decay rate is the Ji,4 (i.e. 
J2), the exchange interaction along the Fe-Fe square di- 
agonal direction. Ji_4 is very large (i.e. way above the 
1 / i?'^-envelop) and then basically goes to zero for further 
distances (i.e. Ji^i^ = Ji iq « 0). This suggest that the 
origin of Ji^ is probably not the Fermi-surface nesting 
or other long-range exchange interactions but rather lo- 
cal superexchange interactions through As-p orbitals. In 
contrast, the nature of the exchange constants along the 
a— and 5— directions are very different. They are long 
range and decay with l/i?"^, much like in bcc-Fe^. 

The bottom panel in Fig. 7 shows the calculated ex- 
change parameters Jij{R) obtained from linear response 
perturbation approach with rigid-spin approximation. 
Similar to results from spin-flip method, the exchange in- 
teraction along the Fe-Fe square is antiferromagnetic and 
very short range while along the a— and 6-directions it 
has oscillatory character with sign change and slow decay. 
The most important difference between the spin-flip and 
perturbation methods is the nearest-neighbor exchange 
interaction along the stripe direction, i.e. Ji,3. While 
spin-flip method gives this interactions as antiferromag- 
netic the linear response theory suggests it is weak but 
ferromagnetic. As discussed above and shown in Fig. 5, 
this difference between two methods is a nice evidence 
that the spin-interactions along the stripe direction is 
probably not a classical Heisenberg. In conclusion, due 
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FIG. 7: (color online) The calculated exchange constants 
( S in JS^ is taken to be 1) along three different 
paths as shown in the inset. Top panel shows the results 
from direct spin-flip method while the bottom panel shows 
results from linear response perturbation method within rigid- 
spin-approximation. Both methods give strong and antiferro- 
magnetic short range interaction along the diagonal direction 
while they differ for the parallel aligned spin-direction (i.e. 
b-axis). The dashed curves and the shaded area in top panel 
shows that the exchange interactions along the a and b axes 
are oscillatory with an envelop decaying as l/B?, much like 
in the bcc-Fe^. 

to slow 1 /R^ decay of Js along the stripe direction and its 
different sign obtained from spin-flip and linear response 
methods strongly suggest that the spin-interactions along 
the stripe direction is unusual, long-range and probably 
related to Fermi-surface nesting. 

We note that our resuhs from perturbation theory are 
in agreement with the earlier calculations reported from 
Pickett's group in Ref,46 where the spin-exchange con- 



stants were calculated in various 122 systems using linear 
response theory and small but ferromagnetic exchange in- 
teractions are found along the parallel-spin direction. We 
emphasize that since the exchange interaction along the 
stripe direction (i.e b-axis) is ferromagnetic, the frustra- 
tion is totally removed and there is nothing special for 
the ratio Jf/2J2 being « 1 anymore. Finally, we note 
that there is an other interesting report by Balashchenko 
and Antropo"\«^ where the exchange interactions are cal- 
culated in real space as a function of As z-position using 
linear response theory. As expected, they observe mag- 
netization changes from 1.3/xb to 0.31^b with a small 
displacement of arsenic z-position which in turn changes 
the Fe-As bond distance by about 0.1 A. However the 
calculated exchange parameters do not agree with those 
reported from Pickett's group even though both groups 
use linear response theory. They found that all the 
major interactions are antiferromagnetic with significant 
anisotropy along the parallel and antiparrallel spin di- 
rections. This anisotropy is found to be very sensitive 
to the As z position. However for all z-positions studied, 
the exchange along the stripe direction is always antifer- 
romagnetic. The authors also conclude that the inter- 
actions are long range in agreement with our results for 
the spin directions along the a— and 6— axis. However 
we emphasize that from both methods used here, we find 
that the diagonal spin interaction (i.e. J2) is very strong 
and antiferromagnetic for only the nearest neighbor and 
then basically becomes zero for further distances. 

Since two methods shown in Fig. 7 give opposite spin- 
interactions along the stripe direction (i.e. parallel spin- 
direction), an interesting picture emerges from these two 
calculations. At high temperature where the system is 
paramagnetic or close to magnetic ordering , we find 
that the major magnetic interactions are antiferromag- 
netic and frustrated. The exchange interactions along 
a- and b-directions are comparable to each other as it 
should be due to tetragonal symmetry of the paramag- 
netic state. As the system orders and the spin-fiip ex- 
citations become more and more spin-oscillations as in 
the case of spin- waves, the band structure is modified 
(i.e. the symmetry between dxz and dyz is broken) yield- 
ing very different exchange interactions along a- and b- 
directions. In the low-temperature limit where the spins 
are pretty much fixed, we expect the perturbation results 
valid and should replace the spin-flip results. 

Due to non-Heisenberg nature of the spin-interactions 
that we obtain here, one has to be careful in modeling 
these systems using a Heisenberg like model. For a given 
spin-configuration, it is probably OK to use a Heisenberg 
model to study low-energy excitations in that configura- 
tion (such as spin- waves in SDW ordered state) . However 
if one wants to study the whole phase diagram as a func- 
tion of temperature, a simple Heisenberg model with a 
fixed Js is clearly not appropriate. As the spins rotate 
or change configurations, one will obtain totally differ- 
ent exchange constants. Hence in this case, one needs to 
go beyond the classical Heisenberg model to treat the Fe 
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d— orbitals and spin degrees of freedom on equal foot- 
ing. Hence a Kugel-Khomskii like Hamiltonian could 
be more appropriate for Fe-pnictide systems. Recently 
a complicated phase diagram of spin-orbital ordering of 
Fe-pnictide has been studied in Ref 49 and we refer the 
reader to this study for details. 

We finish this section with a brief discussion of the 
experimental determination of the exchange parameters 
in Fe-pnictides from inelastic neutron scattering exper- 
iments. Even though, Heisenberg model is not appro- 
priate as discussed above, it is probably good enough 
to describe the spin-wave excitations at low tempera- 
tures where the spins are not flipping and therefore the 
exchange constants are not changing. A minimal spin 
Hamiltonian for Fe-pnictide may be written as 



i a—a^b,d,c 



J a Si ■ S, 



(1) 



where the i-summation runs over all Fe-spins and Ja 
{a = a, b, c, d) are the exchange interactions along the 
antiparallel spin direction a, along the parallel spin- 
direction b, along the c-axis, and finally along the Fe- 
Fe square diagonal, respectively. The last term is the 
easy- axis single ion anisotropy which originates from 
spin-orbit coupling. Usually it is small but here due to 
very large Ja and J^, it is important to keep this term 
which can give large spin-gap at gamma that is propor- 
tional to \/D X {Ja + 2Jd)- The spin-wave spectrum of 
this Hamiltonian can be easily obtained by considering 
the AF2-spin configuration as helimagnetic ordering with 
modulation wavevector Q = (tt, 0, tt) such that the spins 
are ordered antiferromagnetically along the a— and c— 
axis and ferromagnetically along the 6-axis. We note 
that the unit cell of the helimagnetic spin-structure is 
twice smaller than the chemical cell along all directions 
(i.e. aM = aojl^bu = bol'2., and cm = Co/2). Within 
this helimagnetic description of the cell, the above model 
Hamiltonian has a single spin-wave mode 



.;(q) = 2S^IAI-Bl 



Ja - Jb[l - cos(gfc)] +2Jd + Jc + D 



(2) 



Ja coa{qa) + 2Jd cos(ga) cos{qy) + Jc cos(9c) 



and the zero-temperature inelastic structure factor 
(which is proportional to inelastic neutron scattering in- 
tensity) 



5(q,c.) 



Uq-Bq 

A„ + B„ 



5{UJ - LOq) 



(3) 



After having determined the cigen-spectrum of our 
model Hamiltonian, we now discuss how one may mea- 
sure them using inelastic neutron scatting. Because of 
the tetragonal symmetry of the paramagnetic phase, dur- 
ing the SDW transition the single crystal samples will 
probably have orthorhombic twinning which we call ab- 
domains. In one domain the stripe direction (i.e. parallel 
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FIG. 8: (color online) The afe— domain-averaged spin-wave 
spectrum along (q,0,h) and (q,q,h) directions (h=0, 1). The 
black and gray (red) curves are due to spin- waves coming from 
domain-a and domain-b where the stripe direction is along the 
a- and b-axis, respectively. Note that with a proper constant 
momentum and/or energy scans, it should be possible to re- 
solve two modes from two domains and uniquely determine 
the sign of exchange parameter J;,. On the right, we give 
the analytical expressions for the spin-wave gaps at differ- 
ent Brillouin zone boundaries. The plots are obtained using 
the following values (all in meV); Jd = 40.0, Ja = 20, Jc = 
5.0, _D = 0.015. 



aligned spin direction) will be along the a-axis while in 
the other domain it will be along the b-axis. In other 
words, even though we have a single crystal sample, we 
can not distinguish the wave- vector q = {qa,qb,qc) from 
q — {qb,qa,qc)- Hence, from inelastic neutron scatting 
one will therefore probe the superposition of the spin- 
wave modes at these two wave- vectors, i.e. uj{qa, qb, Ic) + 
(-u{qb, qa, Qc)- As we shall see below, due to very different 
ordering of spins along the a— and b— axes, the super- 
position of these two modes can be, in principle, resolved 
at high energies and in particular at zone-boundaries. 

In Fig. 8 we show the a6— domain averaged spin-wave 
spectrum, i.e. uj{qa,qb, qc)+i^{qb, Qa, Qc), along various di- 
rections. For the exchange parameters, we consider three 
cases; isotropic model ( Ja = Jb), anisotropic model (i.e. 
Jb — Ja/2), and finally ferromagnetic model {Jb < 0). 
The values of the exchange parameters used in the cal- 
culations are given in Fig. 8. From Fig. 8, it is clear 
that contributions due to a— domain (red-curve) and b— 
domain (black curve) are quite different and could be 
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FIG. 9: (color online) The contour plot of calculated S{q,ij) 
for ISO, aniso, and ferro-model for Jf, along (q,0,0)-direction. 
The energy scans at constant q=(l/2,0,0) for the three mod- 
els are indicated by dashed-line and shown in the bottom 
right panel. Note that despite the presence of orthorhombic- 
twinning, the spin-wave spectrum of the three models are 
quite distinct suggesting that one may uniquely determine 
the sign of Jb- All these are based on the assumption that 
the Stoner-continuum will not over-damp the spin-wave spec- 
trum. 



resolved experimentally. Interestingly for the isotropic 
model (Jq = Jb), the spin-wave spectrum is the most 
anisotropic near the zone-center and near the middle of 
the (1/2, 0, 0) as shown in Fig. 8. This is due to the fact 
that the spins are aligned antiparallel and parallel along a 
and b axes, respectively while their interactions are forced 
to be equal and antiferromagnetic. This gives very dif- 
ferent dispersion along the a- and 5— axis and therefore 
one can resolve the two peaks in the a6— domain averaged 
spectrum. From Fig. 8, it is clear that near the (1,0,0), 
the isotropic model have two modes near zero energy 
while the anisotropic and ferromagnetic models have one 
mode that has a huge spin-wave gap that is proportional 
to y/ {Ja — Jb){2Jd — Jb)- The analytical expressions for 
the spin gaps at various Brillouin zone boundaries are 
given in Fig. 8. Finally, we note that {qqO) direction is a 
special one where we should see only a single mode which 
has a tiny gap at (100) if the interactions are isotropic and 
a gap about the half of the maximum spin-wave energy 
for the anisotropic model. For the ferromagnetic model, 
the gap is almost the same as the maximum spin-wave 
energy (see left bottom panel in Fig. 8). 

In Fig. 9, we show the a6— domain averaged inelastic- 
structure factor to get an idea about the intensities of 
the modes from two domains. The calculated spectrum 
is convoluted with Gaussian and 5% energy resolution. 
Fortunately the intensities from both domains are com- 
parable and therefore it should be feasible to detect these 
modes in a multi-domain sample. As an example, in 
Fig. 9, we show energy scans at constant wavevector 



Q = (1/2,0,0) for the three modes, which show very 
distinct spectrum for each model. 

As a final note, we point out that in practice, the ob- 
servation of spin-waves in Fe-pnictide could be problem- 
atic at high energies due to strong spin- wave damping by 
Stoner-continuum. However very recent inelastic neutron 
scattering measurements on Cal22 system^^ has revealed 
steeply dispersive and well-defined spin waves up to an 
energy of ~ 100 meV. Unfortunately the resolution and 
the quality of these recent data are not good enough to 
carry out a detailed analysis as discussed here in order 
to determine the sign of Jj*. We hope that in the near 
future there will be more experimental data with better 
resolution and quality. 



IV. STRUCTURAL PHASE TRANSITION IN 
FE-PNICTIDES 

We next discuss the implication of the magnetically 
frustrated AF2 configuration on the structural distortion 
which is shown to be a common feature of 1111 and 112 
parent compoundsi^. Experimentally it has been demon- 
strated that magnetic SDW ordering and the tetragonal- 
orthorhombic distortion are closely coupled-^. Interest- 
ingly for 1111 systems (i.e. LaOFeAs)^^ , the tetragonal 
distortion first takes place about 20-30 K higher in tem- 
perature than the magnetic transition. The transition for 
1111 system is somewhat weak. On the other hand, for 
the 122 systems (such as BaFe2As2), it is clear that both 
transitions occur at the same time. They are more or 
less strong first-order type and one study suggests that 
the structural distortion is proportional to the magnetic 
order parameter (rather than its square)^. 

In order to establish a connection between the struc- 
tural distortion and the magnetic ordering, we calculate 
the total energy as a function of the 7 angle as shown 
in the inset to Fig. 10 for NM, F, AFl, and AF2 spin 
configurations (see Fig. 3). When 7 = 90°, we have the 
original tetragonal cell. Once the 7 deviates from 90°, 
the original x structure (shown as dashed line) 
is no longer tetragonal but orthorhombic (i.e., the cell 
length along a and b axes are no longer equal). We note 
that for 7 = 90°, the orbitals dxz and dyz are degenerate 
and therefore one may think that the system is subject 
to symmetry lowering for reasons similar to those in a 
Jahn- Teller distortion. However as shown in Fig. 10, we 
do not see any distortion for any of the NM, F, and AF 
configurations. 

The total energy versus 7 plot shown in Fig. 10 clearly 
indicates that only AF2 ordering distort the structure 
with 7 = 91.0°, which is in good agreement with the ex- 
perimental value of 90.3°. The FP-LAPW method with 
LDA approximation we get M=0.48 fis which is in ex- 
cellent agreement with the experimental value of 0.35 
fiB- Hence all electron-calculation is able to reproduce 
both the observed structural distortion and the small Fe- 
moment simultaneously. On the other hand, PW calcu- 
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FIG. 10: (color online) The total energy per cell versus the 
angle 7 for non-magnetic (NM), Ferromagnetic (F) and two 
antiferromagnetic (AFI and AF2) spin-configurations. Note 
that only the AF2 spin configuration yields structural distor- 
tion. The inset shows that as 7 increases, the ferromagnetic 
aligned Fe ions (i.e., Fei-Fe2) get closer while the antiferro- 
magnetically aligned ions (i.e., Fei-Fes) move apart, breaking 
the four-fold symmetry and thus the degeneracy of the d^z 
and dyz orbitals. For the AF2, the solid and dashed lines 
are from pseudo potential plane wave (PW) and FP-LAPW 
calculations, respectively. 



lations with GGA approximation (dotted line) give good 
structural parameters when spin-polarization is allowed 
but then the calculated moment is too large. Somehow 
the magnetic ground state is over-stabilized in PW calcu- 
lations. Therefore one has to be careful in studying the 
magneto-elastic couplings by PW methods, which will 
be underestimated. The net energy gain by the struc- 
tural distortion shown in Fig 10 is about 12 meV per 
cell, which is of the same order as the temperature at 
which this phase transition occurs. We also considered 
two types of AF2 where the spins along the short axis are 
aligned parallel or anti-parallel (see Fig. 3c-d). These 
two configuration are no longer equivalent. According 
to our calculations the configuration in which the spins 
are ordered parallel along the short-axis is the ground 
state. This prediction2i has now been confirmed by neu- 
tron scattering measurements for several 122 systems^^, 
giving us confidence that first-principles calculations de- 
scribe many fine details of Fe-pnictide systems accurately. 

Even though full structural optimization discussed 
above clearly shows that AF2 ordering gives rise to ob- 
served orthorhombic distortion with the long axis along 
the anti-parallel spin direction, it is not obvious why this 
happens. Naively one would expect the opposite, i.e. 
parallel spin-direction is the frustrated one and therefore 
it should get longer to decrease the antiferromagnetic in- 
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FIG. 11: (color online) Top: The paramagnetic portion of 
total energy (Ep) versus the tetragonal lattice parameter a 
(a and b are taken to be equal). The inset shows the SDW 
ordering, the relevant exchange constants and the fit to the 
total energy. Bottom: The linear expansion of the exchange 
constants with respect to lattice deformation. The tetragonal 
cell is distorted along a-direction while the b-axis is kept con- 
stant Bit b — ar = 5.703 A. Points are the actual calculations 
and solid lines are the linear fit. 



teractions. This would be the case if we had Cu-O-Cu 
linear bond. However in the Fe-pnictide case, the arsenic 
ions are not directly between two Fe-ions (see Fig. 2) but 
they are above/below the Fe-squares. Hence when Fe-Fe 
distance is increased along a-axis, the Fe-As-Fe bond an- 
gle also increase along this direction while Fe-As distance 
does not change at first order. In fact. As ions are pulled 
down to keep the rigid FeAs bond fix and increase the 
angle. However the net effect of the increase in angle on 
the exchange interaction itself is not easy to predict due 
to complicated nature of the problem. 

To determine if the exchange interaction Jf increases 
while Ji decreases during the orthorhombic transition, 
we calculate the lattice dependence of the exchange pa- 
rameters, i.e. dJ/da, and then develop a simple spin- 
Peierles like model which successfully explains both the 
transition and the lattice parameters of F, AFf and AF2 
spin configurations within a unified picture. 

In order to obtain the lattice-parameter dependence of 
the Js, we calculate the energies of F, AFI and AF2 con- 
figurations for different tetragonal lattice parameters and 
then extract the Ep, the paramagnetic portion of the to- 
tal energy as a function of tetragonal cell a. This is shown 
in the top-panel of Fig. 11. The tetragonal elastic con- 
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stant is then extracted by fitting the total energy to a 
quadratic form as shown in the figure. After having de- 
termined the paramagnetic contribution of the total en- 
ergy, we then apply tetragonal to orthorhombic distortion 
by varying the lattice parameter along a-direction. For 
each distortion, we then calculate the energies of F, AFl 
and AF2 spin-configurations. For each spin-configuration 
the internal atomic coordinates such as As z-position are 
always optimized. Using the energy expressions given in 
Fig. 3, we then extract the Ji and J2 which are given in 
the bottom panel of Fig. 11. We note that Jf increases 
and Ji decreases linearly as the a-axis is elongated. Sim- 
ilarly, the J2 always increase with increasing lattice pa- 
rameters near the paramagnetic optimum tetragonal cell 
parameter ar as shown in Fig. 11. Hence it seems that 
both Ji and J2 increases with increasing bond-angle. 

Now, using the lattice parameters dependence of 
Ji, Ji, and J2 calculated in Fig. 11, one can easily obtain 
the new lattice parameters of the magnetic cell for ferro, 
AFl and AF2 spin configurations which arc summarized 
below: 



F ^a = b 



AFl a = b = 



AF2 a 



25 a- (3 2 ^ 

a-T - ( 1 ) < S > 

Co Co 

5.6771 (5.6251) 

ax — [ ) < o > 

5.6991 (5.7011) 



,25 a + 13 

AT + (— + 

5.739I (5.7341) 



a-T + ( 



25 a + (3 



eo eo 
5.6961 (5.6681) 



)<S^> 



(4) 
(5) 
(6) 
(7) 



Here the numbers given in parentheses are the results 
from self consistent full cell relaxation calculations. From 
above results, it is clear that this simple model explains 
most of the observed features such as F and AFl order- 
ing does not distort the lattice. Ferro magnetic cell has 
the smallest lattice parameter because the system wants 
to make both Ji and J2 smaller as they are antiferro- 
magnetic and we are forcing the spins ferromagnetically 
ordered. Similarly, for AFl ordering, we get competitions 
between Ji which wants to increase the lattice while the 
J2 term is still not satisfied and therefore it wants the 
lattice shrinks, yielding a lattice parameter larger than 
ferro but smaller than AF2. In the case of AF2, the J 2 
is totally happy and want to increase the lattice. Our 
model for AF2 nicely predicts orthorhombic distortion 
with right lattice parameters. 

Finally we note that all of our discussion given above 
is basically based on a single FeAs layer without inter- 
plane interactions. In reality these systems order three 
dimensionally at the structural phase transition. This 
raises a questions; is there any correlation between the 
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FIG. 12: (color online) Top: The structural phase transition 
temperature Tstr versus FeAs interplane distance for M122 
(M=Ca, Eu, Sr, and Ba) and Lallll systems. Bottom: The 
calculated interplane exchange interactions J± (in meV) ver- 
sus FeAs interplane distance. For the Lallll system, the 
energy difference between AF and F magnetic configurations 
(shown in the inset) is too small to get an accurate number 
but it is less than 0.2 meV. Experimentallyi^ J± in 1111 sys- 
tems is also found to be very weak with both positive and 
negative sign depending on the rare-earth in the system. The 
inset shows the ferro and antiferro alignments of the FeAs 
planes in the unit cell. Here Ep and Eaf are the total en- 
ergies per cell (i.e. 8 Fe atoms). Note that while J± drops 
sharply with increasing Fe-Fe interplane distance, Tstr first 
increases and then drops with increasing Fe-Fe distance. 



structural phase transition Tstr and the inter-plane mag- 
netic interaction J± between FeAs planes. Normally the 
larger the J± the higher the Tstr should be. To answer 
this question, we have calculated J± for different 122 sys- 
tems and plot it as a function of Fe-Fe plane distance in 
Fig. 12. From this figure it is clear that structural phase 
transition temperature Tstr does not have a monotone de- 
pendence on Fe-Fe distance while J± does. This suggest 
that there must be an other factor which effects the Tstr 
in an opposite way. One possibility could be the shear 
modulus of the system. If the c-axis is short, then the 
chemical bonding between FeAs planes would be stronger 
making the structural distortion difhcult. On the other 
hand, the smaller c-axis means the larger J± which means 
larger in-plane AF2 correlation that drives the structural 
phase transition. Since the effects of shear-modulus and 
are opposite, T^t^ could have a maximum at a partic- 
ular FeAs interplane distance as indicated by the experi- 
mental data. Currently we are developing a phenomeno- 
logical Landau-theory along these lines to address the 
remaining issues discussed above and our results will be 
published elsewhere^. We note that there have been al- 
ready a large number of theories^^i^i^iS to describe this 



12 



coupled structural and magnetic phase transitions. We 
refer the reader to these studies for details^SiSliSEiS^. 



V. GIANT MAGNETO-ELASTIC COUPLING IN 
FE-PNICTIDES 

In order to get a general understanding giant magneto- 
elastic coupling present in iron-pnictides, here we con- 
sider one example of each class of pnictides; namely 
CaFe2As2 for the 122 system with the smallest Ca-ion 
available and the BaFe2As2 with the largest metal Ba. 
For the 1111 system, we consider LaOFeAs. We also 
study a doped 122-system, i.e. Nao.5Cao.5Fe2As2. Since 
m our %/2 X V2-cell we have four chemical formula, we 
consider a supercell where two Na and two Ca are or- 
dered. For each given system, we have performed full 
structural optimization including the lattice parameters 
and the atomic positions. We consider our optimization 
is converged when the maximum force on each atom is 
less than 0.005 eV/A and the pressure is less than 0.1 
kbar. We have performed the full structural optimiza- 
tion for non-magnetic (NM), i.e. "non-spin polarized", 
checkerboard antiferromagnetic (AFl) and stripe (AF2) 
spin configurations. Our results are summarized in Ta- 
ble 1. As expected, the ground state for all four systems 
is the stripe AF2 phase and the optimized parameters are 
in good agreement with the experimental data at ambient 
conditions. 

The most striking and surprising finding listed in Ta- 
ble 1 is the giant dependence of the optimized c-lattice 
parameter on the spin-configuration considered. For the 
case of CaFe2As2, we note that AFl configuration is the 
next stable state (after the AF2) but the c- value is signifi- 
cantly reduced; 11.63 A versus 10.60 A for AF2 and AFl 
spin configurations, respectively. This difference in c is 
not due to different AFl and AF2 spin configuration but 
due to different Fe-spin state in AFl and in AF2. The 
difference is even larger, when the Fe-magnetism is to- 
tally ignored (i.e. non-spin polarized calculations). The 
optimized z-value for NM-state is 10.39 A, which is 1.34 
A shorter than the experimental value at ambient pres- 
sure. We note that the optimized lattice parameters, 
a=5.65 A and c=10.39A for the NM phase are in rea- 
sonable agreement with the neutron data in the collapse 
phase (a=5.8 A and c=10.6 A)^. Hence from these 
results and the recent experimental observation of the 
coUapse-T phase, one can reach the conclusion that the 
Fe-moment should be present in 122 systems at all times 
at ambient pressure. This is because we know that these 
systems are ordered in an AF2 spin configuration when 
they are not doped or superconducting. When they are 
doped and superconducting we do not see huge changes 
in their c-lattice parameters. This indicates that even 
though the AF2 long range ordering is destroyed with 
doping or we are at temperatures above the magnetic or- 
dering transition Tjv, the Fe-spin should be present in 
the system. Otherwise, we should see the expected huge 



TABLE L Various optimized structural parameters for NM, 
AF2, and AFl spin configurations, respectively. Note that 
the c-axis from non-spin polarized and AFl-spin configura- 
tions are significantly smaller than the experimental data at 
ambient conditions. The zero of energy is taken as the en- 
ergy of the NM-case. The experimental data are taken from 
Refs^Hliill^. *The AFl configuration goes to NM during 
structural optimization for Cao.5Nao.5Fe2As2 . 





a 


b 


c 


As(z) 


dpeAs 


MFe 


|E (meV) 


CaFe2As2 


NM 


5.63 


5.63 


10.39 


0.36251 


2.309 





0.0 


AFl 


5.65 


5.65 


10.60 


0.36440 


2.338 


1.3 


-16 


AF2 


5.61 


5.48 


11.61 


0.36695 


2.367 


2.2 


-100 


Exp. 


5.68 


5.68 


11.73 


0.3665 


2.370 


1.0 










Cao.5 


Nao.5Fe2As2 






NM 


5.59 


5.59 


10.52 


0.36284 


2.31 





0.0 


AFl* 


5.59 


5.59 


10.52 


0.36284 


2.31 





0.0 


AF2 


5.43 


5.53 


12.05 


0.36536 


2.382 


2.4 


-97 


Exp. 


5.42 


5.42 


11.86 






0.0 




BaFe2As2 


NM 


5.58 


5.58 


12.45 


0.3479 


2.319 





0.0 


AFl 


5.64 


5.64 


12.73 


0.35231 


2.382 


2.1 


-80 


AF2 


5.70 


5.59 


12.83 


0.3549 


2.408 


2.4 


-169 


Exp. 


5.52 


5.52 


13.02 


0.3545 


2.397 


1.0 




LaOFeAs 


NM 


5.64 


5.64 


8.59 


0.35944 


2.332 





0.0 


AFl 


5.69 


5.69 


8.71 


0.35128 


2.393 


2.1 


-86 


AF2 


5.67 


5.73 


8.72 


0.34860 


2.409 


2.4 


-190 


Exp. 


5.70 


5.70 


8.737 


0.3479 


2.407 


0.35 





reduction in the c-axis. In other words, the iron-pnictide 
system should be considered as paramagnetic i.e. with 
the on-site non-zero iron moment without any long range 
order. We note that the paramagnetism is different than 
the " non- magnetic" case that we consider in our DFT 
calculations where we force equal up and down spins in 
each orbital. The non-spin polarized calculations should 
not be considered as a model for the paramagnetic sys- 
tem. With the standard density functional theory, there 
is no way to treat a paramagnetic system (i.e. we have 
the spin-degrees of freedom at each site but no long range 
order). 

Fig. 13 shows an energy reaction path without any en- 
ergy barrier for the collapse of CaFe2As2 with the lost 
of Fe-spin. It is surprising that the c-lattice parameter 
reduces from 11.7 A to 10.4 A but yet the total en- 
ergy change is only 0.3 eV per cell (i.e. four Fe ions). 
One wonders how the atoms are rearranged during the 
collapse of c-axis. Do the FeAs planes remain rigid and 
get close to each other? Does this collapse have anything 
to do with the effective radius of Fe-ion for different Fe- 
spin states? Our answer to these questions is simply no. 
First of all, the collapse of the c-axis happens rather uni- 
formly. There is significant and comparable decrease in 
the height of the Fe-As and As-Ca-As planes, indicating 
that the whole lattice almost uniformly shrinks. In our 
earlier study^^ we traced down the origin of this huge c- 
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Collapsed Tetragonal Phase (cT) 




0123456789 10 
Reaction Path ( T-phase — > cT-phase) 



FIG. 13: (color online) Total energy along a reaction path, 
showing that the Cal22 tetragonal phase goes to coUapsed- 
tetragonal phase without any energy barrier during non-spin 
polarized structural optimization. The energy gain due to 
1.34 A c-axis coUapse is about 0.3 eV. The insets show the 
initial (T-Phase) and final (cT-phase) Cal22 structures with 
relevant bond-distances (in A) and angles (in degrees). Note 
that in the cT-phase the change in the height of the FeAs 
(ALpeAs) and the As-Ca-As (ALAi,As)planes are comparable, 
indicating a uniform compression of the whole lattice. 



axis collapse to large As- As interaction between adjacent 
FeAs planes. 

In order to demonstrate that there are large hybridiza- 
tion between As ions in the Cal22 system, we show the 
contour plots of the relevant orbitals in Fig. 14. Again 
the contour plots in the T- and cT-phases are quite simi- 
lar. It is very clear that the As ion below the top Fe-plane 
makes a bond (or hybridizes) with the arsenic ion which 
is above the lower Fe-plane. Hence this overlap of the 
As- As along the c-axis makes this system quite isotropic 
and far from being layered system. From Fig. 14 it is 
clear that the As- As bonding along the c-axis got signifi- 
cantly stronger. According to bond-population analysis, 
the bond strength increased almost twice. Due to close 
proximity of the As ions in adjacent Fe-layers, the ob- 
servation of the As- As interaction is probably not that 
surprising. What is surprising is to see that there is al- 
most the same type of hybridization between two arsenic 
ions on the same Fe-plane as shown in the bottom panel 
of Fig. 14. 

Since we have shown that the As ion above the Fe- 
plane has a strong overlap with the As ion below the same 
iron plane, their interaction is automatically increased as 
the Fe-As interaction decreases due to decrease in the 
Fe-moment which changes the chemistry of the Fe ion. 
Therefore, we have now a mechanism which explains why 
the As ion z-values get shorter with the decreasing Fe- 
moment. Our mechanism also explains why we see a 




FIG. 14: (color online) Contour plot of one of the orbitals 
which is responsible for the discovered As- As covalent bonding 
for T-phase (top-left) and cT-phase (top-right), respectively. 
Note that the As- As bonding present in both phases is much 
more significant in the cT-phase. The bottom panel shows 
an other orbital on a slice along (110) plane, indicating clear 
hybridization between intra-As atoms below and above the 
Fe-plane in the T-phase. 

smaller reduction in the c-axis for the LaOFeAs than the 
122 system as listed in Table 1. The reduction in the 
c-axis in the LaOFeAs system is due to the intra-plane 
As-As interaction only since there are no two adjacent 
FeAs planes to interact with each other as in the case 
of Cal22. Our theory also predicts that for larger ions 
like Ba, we should see less c-reduction because the As-As 
distance between two adjacent planes are now larger due 
to larger ionic radius of Ba. In Table 1, we also show 
that similar c-reduction with Fe-spin occurs in the doped 
Nao.5Cao.5Fe2As2 system as well. 

Since our results suggest that Fe-magnetism is totally 
lost in the cT-phase, one wonders if the ~ 12 K super- 
conductivity observed in some experiments in the vicin- 
ity of the collapse cT-phase of CaFe2As2^ can be ex- 
plained by conventional electron-phonon (e-ph) coupling? 
In order to address this question, we have calculated 
phonon spectrum and Eliashberg function from linear 
response theory^^. We used basically the same method 
and the equivalent parameters that are used in ReflBVl 
for LaOFeAs. Our results are summarized in Fig. 15 
and very similar to those for LaOFeAs. We obtained 
a value of electron-phonon coupling A = 0.23 and the 
logarithmically average frequency uiog = 218 K, which 
gives Tc = 0.6 K using the Allen-Dynes formula with 
/i* = (i.e. an upper bound for Tc). Hence, if the 
~ 12 K superconductivity observed in some experiments 
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00 (cm ) 

FIG. 15: (color online) Phonon density of states (DOS), 
Eliashberg function {a'^F{Lj)) and the frequency-dependent e- 
ph coupling \{uj) (dashed line) for CaFe2As2 in the cT-phase. 



can be confirmed to be actual bulk superconductivity in 
the cT-phase, it would mean that the mechanism of su- 
perconductivity in the cT-phase of CaFe2As2 is likely un- 
conventional and it has nothing to do with Fe-magnetism 
which is not present in the cT-phase. On the other hand, 
if the future experiments totally rule out that the cT- 
phase of Fe2As2 does not show bulk superconductivity 
then it would mean that the Fe-magnetism is required 
for the superconductivity in these systems. Because of 
these reasons, the cT and T-phases of CaFe2As2 provide 
us an invaluable opportunity to study the same system 
with and without Fe-magnetism. We hope that there 
will be more focus on the superconducting properties of 
CaFe2As2 system under pressure to resolve the outstand- 
ing issues about the ~ 12 K superconductivity observed 
in some experiments. 

In conclusion, we have revealed surprisingly strong As- 
As interactions for both intra- and inter-plane arsenic 
ions. The strength of this interaction is controlled by 
the Fe-As chemical bonding. Reducing the Fe- moment, 
reduces the Fe-As bonding, which in turn increases the 
As- As interaction along the z-axis, causing arsenic atoms 
on opposite sides of Fe-square lattice to move towards 
each other. This explains the high sensitivity of the z- 
atom positions and the large reduction of the c-axis with 
the Fe- magnetic moment. From visualization of the elec- 
tronic orbitals, we show that the 122 systems should not 
be considered as layered systems since the As- As inter- 
plane interaction is as strong as the intra-plane As- As in- 
teraction. The 122 system are in fact quite 3D isotropic 
than one initially thinks. Since at ambient pressure, we 
do not observe large c-axis drops in the superconduct- 
ing samples, we conclude that the Fe-magnetic moment 
should be present at all times in these systems, at least in 
122 materials such as CaFe2As2. In other words, the iron- 
pnictide system should be considered as paramagnetic 
(i.e. Fe-moment is present without long range order). 
Non-magnetic treatment of Fe ions changes the chem- 
istry significantly and is not suitable for description of 
these systems at ambient pressure. The giant coupling of 
the on-site Fe-magnetic moment with the As- As bonding 
that we have discovered here may provide a mechanism 



for the superconductivity. Since earlier electron-phonon 
(ep) coupling calculations^ were done ignoring the Fe- 
moment, our results raise some questions about the va- 
lidity of these calculations. Currently we are extending 
such e-ph coupling calculations to magnetic systems us- 
ing finite-displacement method in which the magnetic re- 
sponse of the system to the atomic motion is treated fully 
unlike the standard linear-response perturbation theory. 



VI. MAGNETIC PHONONS 

Whenever a new superconductor is discovered, the first 
thing to do is to check whether the conventional electron- 
phonon (e-ph) coupling can explain the observed tran- 
sition temperature or not. This was also the case for 
Fe-pnictide superconductors. Early electron-phonon cou- 
pling calculations^ based on standard non-spin polar- 
ized perturbation theory indicate that conventional e-ph 
can not explain the observed high temperature in Fe- 
pnictide superconductors. The phonon spectrum and its 
temperature dependence of various 1111 and 122 systems 
have been also extensively studied by inelastic neutron 
scatterin g^^i^^i^° , inelastic xray scatterin g^^i^^i^'^ and by 
nuclear-resonance spectrum which is Fe-specifio2i. These 
studies did not find any significant changes in the ob- 
served phonon spectrum with superconductivity. How- 
ever some of these measurements indicate features which 
are not produced in the standard linear-response non- 
magnetic phonon calculations. For example, Fukuda at 
alJ^. found that the calculated phonon DOS agrees with 
the experimental spectrum provided that the computed 
FeAs force constant is reduced by 30%. Similarly, Reznik 
et aS^ recently observed that in BaFe2As2 system, the 
Ag As mode is around 20-22 meV while there is no such 
feature in the calculated spectrum. Similar observations 
have been made by inelastic neutron scattering^ where 
experimetnal DOS has a nice sharp peak around 20 meV 
while in the calculated spectrum there is nothing in that 
energy range. In addition to these observations, there 
is also a recent Fe-isotope measurements where an iso- 
tope coefficient of 0.4 is observed--. Some anomalous 
electron-phonon interaction in doped LaOFeAs has been 
also reported from first-principles calculation*^. All 
these studies suggest that it is probably too early to rule 
out a possible mechanism based on phonon-mediated su- 
perconductivity in Fe-pnictide systems. 

So far we have shown that the spin-polarized cal- 
culations recover from the failure of non-magnetic cal- 
culations in terms of lattice parameters and internal 
atomic coordinates. Here we show that magnetic cal- 
culations also resolve the most of the outstanding issues 
with the observed phonon modes discussed above. Our 
phonon calculations are done using the plane wave code 
pwscf within finite displacement technique as described in 
RefFrzl. The advantage of direct finite displacement tech- 
nique over the standard linear response theory is twofold. 
The first advantage is that we can do phonon calcula- 
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TABLE II: The symmetries and energies (in meV) of the optical phonons of LaOFeAS in P4/nmm and Cmma phases. The 
energies of the IR-active modes are taken from Ref.T^. The * indicates a significant disagreement. The details of calculations 
and the animations of the modes can be found at http://www.ncnr.nist.gov/staff/taner/laofeas 



r{P4/nmm) = 2 Aig (R) + 4 A2„(IR) + 4 E„ (IR) + 4 Eg(R) + 2 Bi^ (R) 
T{Cmma) = 2 Ag (R) + 2 Big (R) + 4 Bi„(IR) + 4 B2g(R) + 4 B2„ (IR) + 4 B39(R) + 4 B3„ (IR) 


P4/nmm Cmma IR 


P4/nmm Cmma IR 


P4/nmm Cmma IR 


P4/nmm Cmma IR 


P4/nmm Cmma IR 


E„ 7.3 7.4-7.5 - 
Aig 24.9 25.1 
Eg 35.6 35.9- 36.1 - 


A2u 9.9 10.1 12.1 

Big 26.6 26.9 

E„ 34.3 34.6- 34.7 42.0* 


Eg 14.0 14.1-14.2 - 
A2„ 31.2 31.6 30.9 
A2„ 48.6 49.1 53.8 


Eg 17.6 17.7-17.8 - 
E„ 33.7 34.0-34.1 33.2 
Eg 51.6 51.8-52.6 - 


Aig 22.1 22.3 - 
Big 35.2 35.6 - 
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FIG. 16: (color online) The IR data from LaOFEAs sample 
at different temperatures (adapted from Ref.^^), indicating 
very strong temperature dependence for the mode near 340 
cm~^, which is calculated . Our calculations do not produce 
any zone-center phonon near this energy (see Table. 2). 



tions with different atomic displacements (usually, 0.01, 
0.02, and 0.04 A) and determine if there is any anhar- 
monic phonons. For harmonic phonons, the results do 
not depend on the magnitude of the displacement. The 
second and probably the most important advantage is 
that we treat the magnetism and phonon displacements 
equally and self consistently. We have already seen that 
Fe-magnetic moment is very sensitive to the As-z position 
and therefore it is not a good approximation to assume 
that the spins are fixed as the atoms move according to 
a given phonon mode. In our approach, this direct and 
strong interplay of magnetism and structure is treated 
self consistently. 

We first discuss the phonon spectrum of 1111 system 
(i.e. LaOFeAs) in the high (i.e. paramagnetic) and 
low temperature phases and make some comparison with 
available IR-data which is shown in Fig. 16. Since in 
the low-T phase, we have both magnetic ordering and 
tetragonal-orthorhombic structural phase transition, it 
is instructive to study the effect of magnetic ordering 
and structural distortion separately. Hence, we first ig- 
nore the Fe-magnetism and consider phonon spectrum 
when the system is non-magnetic with tetragonal and 



orthorhombic lattice parameters. Our results are sum- 
marized in Table 2. From the the symmetry decom- 
position of the optical phonons in both P4/nmm and 
Cmma phases, we note that the distortion does not in- 
troduce any new IR active modes but rather just splits 
the doubly-degenerate modes into non-degenerate ones. 
However the splitting is quite small; the largest is around 
0.2 meV. This explains why no new modes appear in the 
optical measurements after the transition. We also note 
that the agreement for the energies of the zone center 
phonons with IR data is not as good as one expects. 
In particular, the Eu mode observed at 42 meV is cal- 
culated to be 35 meV, a significantly lower value. In- 
terestingly, this particular mode has a strong tempera- 
ture dependences^ as shown in Fig. 16. This raises some 
red fiags about the possibility of anharmonic phonons in 
these systems. 

Since the As-z position and the Fe-magnetic moment 
are tightly coupled, one may naively expect that the c- 
polarized As A^ mode could be quite anharmonic if we 
perform magnetic calculations. We have checked this first 
by performing a frozen-phonon type calculation where we 
calculate the total energy as the As-atoms are displaced 
along c-axis in LaOFeAs system (see Fig. 17). Please 
note that this is not true Ag mode where both As and La 
atoms are allowed to move along c-axis. We will discuss 
the full phonon calculations next. Figure 17 shows that 
the total energy versus the As-displacement can be fit to 
quadratic equation quite well, suggesting this phonon is 
harmonic. What is interesting is that the non-spin po- 
larized and spin-polarized calculations give very different 
harmonic force constants. Including Fe-spin in the calcu- 
lations softens the force constant by about 20%, consis- 
tent with Fukuda's observation^^ that one has to soften 
the FeAs-force constant to get better agreement with the 
data. However as we shall see below, the renormalization 
of the force constants with Fe-spin polarization is quite 
complicated; almost all of the force constants involving 
Fe and As are basically renormalized. 

We next check if there are other modes in LaOFeAs 
system that could be anharmonic. We have carried out 
full spin polarized phonon calculations (gamma point of 
the \/2 X \/2 cell) with atomic displacements of 0.02 A 
and 0.04 A, respectively. The mode energies are plotted 
in Fig. 18 for both displacements. It is clear that two 
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FIG. 17; (color online) Fe-moment and total Energy as As 
ions are translated along the c-axis as shown in the inset. 
Note that -l-(z-zO) corresponds to As ions moving towards the 
Fe-plane (z0=0.36). Bottom panel indicates that the frozen- 
phonon potential is harmonic and softens about 20% with 
spin-polarization. 
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FIG. 18: (color online) Phonon mode energies versus mode 
numbers from finite-displacement phonon calculations with 
two different displacements, indicating that most of the 
phonons are harmonic. There are some modes with very small 
anharmonicity around 250 cm~^ which correspond to Fe-Fe 
stretching modes. The As c-polarized modes with Ag sym- 
metry are indicated. These modes are harmonic, consisting 
with frozen-phonon energy plot shown in Fig. 17. 



As Ag modes are harmonic. There is some anharmonic- 
ity around 250-300 cm~^ corresponding to in-plane Fe-Fe 
stretching modes. However this type of mode energy de- 
pendence on the atomic displacement is typical and does 
not indicate unusual anharmonicity. Hence, we conclude 
that in Fe-Pnictide system the phonon spectrum is har- 
monic. This is quite opposite to what we have in MgB2 
superconductor-^^ where the in-plane B-B stretching 
mode was found very anharmonic and responsible for the 
most of the e-ph coupling2i. 

We now discuss the effect of the Fe-spin state on 
phonon modes in LaOFeAs. Our results are also very 
similar to 122 systems (in fact the effect of the Fe-spin 
is more pronounced in 122 systems than in 1111 systems 
due to strong inter-plane As- As interactions as discussed 
above). We have carried out three different phonon cal- 
culations. In the first two, we ignore Fe-magnctism and 
just consider tetragonal (i.e. P4/nmm) and primitive 
cell of the orthorhombic space group Cmma (i.e. P2/c). 
The third calculation is fully spin-polarized with the or- 
thorhombic cell parameters. We note that when iron 
spins are considered, the space group of the SDW or- 
dered system is no longer Cmma as reported in neutron 
experiments^*. Cmma is the space group of the non- 
magnetic orthorhombic cell. In our phonon calculations 
this is a very important point since we use only symmetry 
independent displacements to construct the dynamical 
matrix. Using wrong space group, such as Cmma, could 
average out the anisotropic force constants due to SDW 
ordering and give wrong results. We determine that the 
true space group of LaOFeAs when iron spin is consid- 
ered is Pbmb (spg. number=49, origin 1, a-cb). We note 
that even this space group is an approximate since we 
assume that FeAs planes are ordered fcrromagnctically 
along c-axis while they are actually ordered antiferromag- 
netically. However we checked that the energy difference 
between ferro and antiferro ordering of the FeAs-planes is 
too small to have any significant effect on the calculated 
phonon spectrum (see J± in Fig. 12 for Lai 11). Our re- 
sults from these three calculations are shown in Fig. 19 
and Fig. 20. 

Fig. 19 shows how the Arsenic c-polarized A^ mode 
is effected by the structural and magnetic ordering. We 
note that due the LaO-plane in 1111 system, we have two 
Ag modes. In the first one As- and La-atoms move along 
c-axis in phase. In the second one As and La atoms move 
opposite (i.e. out-of-phase) along the c-axis. Hence in the 
out-of-phase Ag mode, the As-La distance changes as the 
atoms move. This causes the out-of-phase Ag mode to 
have slightly higher energy than in-phase mode. We also 
note that As atoms move twice more than La atoms in 
the in-phase mode. This is reversed in the out-of-phase 
mode. Because of this difference, the effect of the Fe- 
spin state is small on the out-of-phase mode. The in- 
phase Ag mode energy is soften by about 10% from 22.3 
meV to 20.18 meV with the iron spin. This is significant 
and indicates large magneto-elastic interactions in these 
systems as we have already seen in the case of CaFe2As2. 
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As Ag (in-phase) 

P4/nmm Non-mag. 
Ag =22.1 meV 

P2/c Non-magnetic 
Ag = 22.3 meV 

Pbmb SDW (spol) 
Ag= 20.18 meV 




As Ag (out-of-phase) 

P4/nmm Non-mag. 
Ag=24.9 meV 

P2/c Non-magnetic 
Ag= 25.1 meV 

Pbmb SDW (spol) 
Ag 24.46 meV 



FIG. 19: (color online) Top panel shows two Arsenic c- 
polarized Ag modes, which are in-phase and out-of-phase 
with respect to As and La motions along c-axis. The bot- 
tom panel shows the mode energies for non-magnetic tetrago- 
nal (P4/nmm), non-magnetic orthorhombic distorted lattices 
(P2/c), and SDW magnetic configuration (Pbmb). 



The effect of the spin-polarization on the whole phonon 
spectrum is shown Fig. 20. The phonon density of states 
are obtained by calculating the dynamical matrix in a 
2-\/2 X 2-\/2 supercell with and without spin-polarization. 
The effect of the Fe-spin on the phonon DOS is signifi- 
cant. As we have already seen, the first effect is softening 
the in-phase As Ag c-polarized mode from 22 meV to 20 
meV. The 2nd largest effect occurs near the 36 meV en- 
ergy range. In the non-spin polarized calculations, we ob- 
tain very strong and sharp feature near 36 meV, inconsis- 
tent with the experimental data. This sharp feature near 
35-36 meV is combination of several modes which are 
both c-polarized and in-plane oxygen and Fe-Fe modes. 
When we have the Fe-spin included in the calculations, 
we soften those modes that involve Fe-Fe stretching by 
about %10, bringing the Fe-modes down to 32 meV where 
the experimental features are observed. 

We also checked the real-space force constants ob- 
tained from both spin-polarized and non-magnetic cal- 
culations. While La and O onsite force constants are not 
effected with Fe-spin, the onsite Fe and As force constants 
are renormalized by about 10-20%. For example, the 
non-spin polarized case gives force constants (in eV/A^) 
(11.06,11.06,8.7) for Fe and (10.48, 10.50, 9.33) for As, 
respectively. The spin-polarized calculations renormalize 
these force constants to (10.8, 8.5, 8.5) for Fe and (9.2,8.8, 
8.6) for As. The change in the force constants are sig- 
nificant. We also observe similar softening up to 20% 
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FIG. 20: (color online) Phonon DOS at optimized orthorhom- 
bic cell (top) and room-temperature tetragonal cell (bottom) 
for non-spin polarized (black) and spin-polarized (red) cases. 
Note that for non-spol case, we have an intense peak at around 
34-35 meV. This intense peak is due to oxygen c-polarized 
phonons (35.50 meV) and oxygen in-plane phonons (34.76 
meV) as well as Fe-Fe in-plane stretching mode (34.88 meV) 
and a mixed As c-polarized mode. When spin-polarization is 
included all modes associated with Fe-Fe stretching and As 
c-modes soften and go down to lower energies. 



in the Fe-As and Fe-Fe force constants. Hence the 10% 
phonon softening of the As Ag mode and Fe-Fe stretching 
modes near 35 meV is due to complicated renormaliza- 
tion of the force constants rather than a simple rescaling 
of a single force constant as suggested by Futuda et alM. 
The important point is that Fe-spin is needed to get the 
observed spectrum even though the measurements are 
done on samples at temperatures well above the Tjv. 

Fig. 20 also shows that using room temperature or op- 
timized lattice parameters gives only slightly different 
phonon spectrum. This is consistent with the observed 
weak temperature dependence of the neutron or inelastic 
x-ray data. Of course, in reality there is a huge differ- 
ence between the room temperature and low temperature 
calculations if we ignore the Fe-magnetism at room tem- 
perature and consider it at low temperature. 

We will finish this section by presenting phonon spec- 
trum of Bao.5Ko.5Fe2As2 for which the difference between 
non-magnetic and magnetic DOS is huge as expected 
from our previous work on CaFe2As2 due to strong As- 
As interactions. In Fig. 21 we compare our calculated 



phonon DOS with the neutron data of Zbiri et al^ on 
BaFe2As2 sample. The authors indicated that a peak 
near 20 meV can not be produced from non-spin polar- 
ized calculations. Zbiri's neutron data shown in Fig. 21 
is fully consistent with the recent inelastic X-ray study 
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FIG. 21: (color online) Generalized Phonon DOS mea- 
sured by inelastic neutron scattering at room temperature for 
BaFe2As2 (bottom black curve)^* and the calculated GDOS 
with (middle blue) and without (top red) Fe-magnetism. Note 
that including Fe-spin in the calculations softens the Fe-Fe in- 
plane and As c-phonons by about %10 and %23 and results a 
DOS that is in perfect agreement with the room temperature 
data (i.e weU above Tn = 140 K). 



of Reznik et al7^ . They have also measure c-polarized 
As dispersion mode near 20-22 meV but their non-spin 
polarized calculations did not produce any peak in the 
energy range of 20-26 meV— . Zbiri et al. called this 
energy range as "pseudo-gap" in the phonon spectrum 
since in the dispersion curves there are no modes in this 
energy range that can give a sharp peak in the DOS. In 
Fig. 21 we show that considering magnetic phonons again 
solves all the mysterious! Our non-magnetic phonon cal- 
culations are in good agreement with Reznik phonon 
dispersion curve^^ as well as Zbiri's results^^ and there- 
fore it can not explain the observed GDOS of BaFe2As2. 
However our magnetic phonon calculations soften the c- 
polarized As mode by 23% and the Fe-Fe modes by about 
%10-14, respectively. Hence the magnetic phonon DOS 
is now in perfect agreement with the experimental data. 
This further supports our theory that Fe-magnetism al- 
ways present in Fe-pnictide systems even at temperatures 
well above T^r . The iron magnetism could be in the form 
of fluctuating SDW type small magnetic domains or it 
could be at the atomic limit of paramagnetic Fe ions (i.e. 
para-magnon) . More study is needed to have a better un- 
derstanding of the Fe-magnetism in these systems. From 
the results presented here it is clear that iron-magnetism 
in the Fe-pnictides is the key factor that controls atomic 
positions, lattice parameters, structural phase transition, 
phonon energies, and most probably the superconducting 
properties as well. 



VII. CONCLUSIONS 

Our conclusions can be summarized as follows: 

• Accurate all-electron fix-spin total energy calcula- 
tions indicate that the ferromagnetic and checker- 
board antiferromagnetic ordering in Fe-pnictides 
were not stable and the stripe Fe-spin configura- 
tion (i.e. SDW) is the only stable ground state. 
Mapping the energies of ferromagnetic, checker- 
board AF and SDW spin configurations on to an 
approximate Jf - Jf- J2 model indicates the presence 
of competing strong antiferromagnetic exchange in- 
teractions in these systems. This suggests that 
magnetism and superconductivity in doped Fe- 
pnictides may be strongly coupled, much like in the 
high-Tc cuprates. 

• The magnetic stripe SDW phase breaks the tetrag- 
onal symmetry, removes the frustration, and causes 
a structural distortion. A simple model based on 
spin-Peierls like transition (i.e. the lattice parame- 
ters dependence of the exchange constants Ji and 
J2) is shown to give the correct amount of lat- 
tice distortion as well as predicts the fine details 
such as in SDW ordering parallel aligned spin- 
direction gets shorter while the anti-parallel aligned 
spin direction gets longer. We show that the struc- 
tural transition temperature, T^tr does not corre- 
late with FeAs-inter-plane magnetic interaction Jx . 
The coupling of the magnetic fluctuations to the in- 
plane strains e^x — ^yy and occupations of orbitals 
dxz — dyz may change the nature of magnetic tran- 
sition to first order and splits the magnetic and 
structural ordering temperatures depending on the 
details of the system^i^i^i^iS. 

• The magnetic exchange interactions Ji_j{R) are cal- 
culated as a function of Fe-Fe distance R from two 
totally different approaches. Both methods indi- 
cate that the exchange interactions along the Fe-Fe 
square-diagonal spin-direction are short range and 
antiferromagnetic. Hence it is tempting to conclude 
that the main diagonal interaction, J2 is superex- 
change type and important contributer towards the 
stabilization of SDW ordering. Along the parallel 
and antiparallel aligned spin-directions, however, 
the exchange interactions have oscillatory charac- 
ter with an envelop decaying as just like in 
bcc-Fe. The major difference between two methods 
is that the first nearest-neighbor exchange interac- 
tion along the parallel spin direction is found to 
be antiferromagnetic in spin-flip method and ferro- 
magnetic in linear-response theory. Assuming both 
methods are accurate, this implies that a simple 
Heisenberg model is not appropriate for the Fe- 
pnictide systems (see Fig. 5). However for a given 
orbital order and spin-configuration, one can still 



19 



use it to describe the low-energy excitations such 
as spin- waves. 

Since the spin-flip method is more appropriate at 
high temperatures above the magnetic phase tran- 
sition, we obtain results that are close to param- 
agnetic tetragonal symmetry (i.e. ~ J^) and 
the system is fully frustrated. As the system or- 
ders, the occupancy of the dxz start differ from 
the occupancy of dyz orbital, giving rise to orbital- 
dependent exchange spin-interactions. At low tem- 
peratures where spin- wave approximation is valid, 
the exchange interaction along the stripe direction 
becomes weak and ferromagnetic. At this point, 
the spin-frustration is totally removed (i.e. all 
magnetic bonds are satisfied). Hence, it would be 
very interesting to determine the exchange interac- 
tions and their temperature dependence by inelas- 
tic neutron scattering measurements. We gave a 
brief discussion about how the spin- wave spectrum 
would appear from a single crystal sample with 
orthorhombic twinning. Due to large anisotropy 
of the spin-wave velocities within the ab-plane, it 
may be possible to resolve the spin waves along the 
a and b directions and determine the sign of the 
Ji exchange interaction despite the orthorhombic 
twinning. 

• We unravel surprisingly strong interactions be- 
tween arsenic ions in Fe-pnictides, the strength of 
which is controlled by the Fe-spin state in an un- 
precedented way. Reducing the Fe-magnetic mo- 
ment, weakens the Fe-As bonding, and in turn, in- 
creases As-As interactions, causing a giant reduc- 
tion in the c-axis. For CaFe2As2 system, this re- 
duction of c-axis with the loss of the Fe-moment is 
as large as 1.4 A, an unheard of giant coupling of 
local spin-state of an ion to its lattice. Since the 
calculated large c-reduction has been recently ob- 
served only under high-pressure, our results suggest 
that the iron magnetic moment should be present 
in Fe-pnictides at all times at ambient pressure. 

• Finally, we showed that Fe-magnetism is also the 
key in understanding the phonons in Fe-pnictides. 
Our magnetic phonon calculations clearly indicate 
that the observed phonon-DOS at room temper- 
ature is much closer to the calculated magnetic 
phonon-DOS rather than non-magnetic one. We 
find that the in-plane Fe-Fe and c-polarized As 
phonon modes are soften by about %23 and %10 for 
122 and 1111 systems, respectively, explaining the 
observed inealastic x-ray data by Fukuda et alJ^ 
and by Reznik et alJ^, and the INS room temper- 
ature GDOS data on BaFe2As2 by Zbiri et 
This finding further supports our theory that the 
Fe-magnetism must present in these systems all the 
time. 

• The main conclusion of our work is that there is a 



giant magneto-elastic coupling in Fe-pnictides. We 
can successfully predict lattice parameters, atomic 
positions and phonons in these systems from first 
principles provided that we always consider Fe-spin 
in our calculations. Since the current electron- 
phonon calculations were carried out without the 
Fe-spin, it is probably too early to rule out el- 
phonon coupling as a possible mechanism. It is very 
important that electron-phonon coupling is calcu- 
lated self-consistently with the Fe-spin and with- 
out the rigid-spin approximation since Fe-moment 
is very sensitive to arsenic motion. Currently we 
are carrying out such calculations and the results 
will be presented elsewhere. 
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APPENDIX A: EXCHANGE INTERACTIONS 
FROM DIRECT SPIN-FLIP METHOD IN A 
LARGE SUPERCELL 

We developed a systematic approach where the ex- 
change parameter between spin-i and spin-j is obtained 
from the total energies of a reference magnetic configu- 
ration and those configurations obtained by flipping the 
spins i and j one at a time and simultaneous flipping of 
both spins. From these four energies, it is possible to 
obtain the exchange constant between spin i and j. We 
note that here we are interested in the isotropic exchange 
interactions. We also do not consider spin-orbit interac- 
tions in our calculations. Hence all calculations are done 
for coUinear spin-configurations. 

In order to extract superexchange interactions up to 
a large cutoff distance, we calculated the total energy 
for various periodic spin configurations based on SDW 
alignment of the z-components of spin {Sz = ±1) with 
a 3-\/(2) X 3\/2 supercell of the LaOFeAs which contains 
144 atoms (36 of which are Fe ions). Since the spin con- 
figuration is the same from one supercell to the next one 
we may write the total energy Ei as 

El = i;o + i^^ J(0,A;;R,0^fc(0)5,(R) 

R k,l 

k,l 

where S(n, R) = Sn is the spin of the nth ion in the 
supercell at R and because of periodicity K{k, I) — 
K{l,k) = J(fc, 0; R). It is obvious that we can 
only expect to determine K{k, I) and not the individual 
J's. However, since the supercell is reasonably large, we 
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can identify the K's with the J at the minimal separa- 
tion. It is also obvious that we can only hope to deter- 
mine K{i,j) for i j, since the energy involving K{i,i) 
depends on (Si)^ = 1 since Si — ±1 for Ni spins. To de- 
termine K{i,j) for i j we calculate four total energies, 
El and the other three corresponding energies when we 
independently change the sign of Si and Sj. When we 
change the sign of Si we get 

E2 = Eo + ^Y. - ^kk]Si[l - 2<5.m](A2) 

where dn,m = 1 if n = m and is zero otherwise. Likewise 
when we change the sign of Sj we get 

E3 = Eo + IJ2 O-^fcil - 2^,, fc]5/[l - 2^j- z](A3) 

and when we change the sign of both spins i and j we 
get 

k,l 

x[l-26,,i][l-25,,i] . (A4) 



Then we construct the quantity X = Ei — E2 ~ E3 + E4, 
to get 

k,l 

[2Sj^k + 2Sjj - 'iSj^kSj.j] (A5) 

Since we require that i ^ j, this gives 

X = AK{i,j)S,Sj , (A6) 

from which we can extract the value of K{i,j). If the su- 
percell is large enough, one can keep only the interactions 
between the nearest neighboring supercell images of the 
spins and therefore the calculated exchange parameter 
can be attributed to spin-interaction between the closest 
pairs of spins of types i and j. We also note that there 
are cases where spin j is at the mid-point between spin 
Si{0) and one of its images at S'i(R). In that case, the 
calculated superexchange constant is twice of the J^.j. 
Similarly there are cases where the spin j is at a point 
where it interacts equally with four images of the spin i. 
In those cases, the calculated J is four times Jij. 
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